Uncovering hidden network architecture from spiking activities using an exact statistical input-output relation of neurons

Identifying network architecture from observed neural activities is crucial in neuroscience studies. A key requirement is knowledge of the statistical input-output relation of single neurons in vivo. By utilizing an exact analytical solution of the spike-timing for leaky integrate-and-fire neurons under noisy inputs balanced near the threshold, we construct a framework that links synaptic type, strength, and spiking nonlinearity with the statistics of neuronal population activity. The framework explains structured pairwise and higher-order interactions of neurons receiving common inputs under different architectures. We compared the theoretical predictions with the activity of monkey and mouse V1 neurons and found that excitatory inputs given to pairs explained the observed sparse activity characterized by strong negative triple-wise interactions, thereby ruling out the alternative explanation by shared inhibition. Moreover, we showed that the strong interactions are a signature of excitatory rather than inhibitory inputs whenever the spontaneous rate is low. We present a guide map of neural interactions that help researchers to specify the hidden neuronal motifs underlying observed interactions found in empirical data.

probabilities. Here, we performed simulations to examine the discrepancies between the theory and simulation caused by these factors.
We traced the spontaneous spiking probability (F 0 ( )) of the LIF postsynaptic neuron by using MATLAB code when the neuron received 5 Hz Poissonian signaling input on top of noisy Gaussian inputs (mean and variance) near the threshold voltage. We used a 5 ms time window to represent the spiking activity and counted the number of bins in which the postsynaptic neuron generated a spike, but the signal did not occur. We divided it by the total number of bins lacking the signal and calculated the spiking probability. The neuron received signaling input at 5 Hz on top of noisy Gaussian inputs, balanced near the voltage near the threshold regime. The membrane time constant was ⌧ m = 10 ms, and the threshold voltage was V ✓ = 5 mV. The simulations ran for duration 10 8 steps with step size 0.1 ms. We repeated the simulation 500 times. Supplementary (e, f) Triple-wise interactions as a function of the amplitude of common inputs in inhibitory-to-trio (e) and inhibitory-to-pairs (f) motifs (filled circles, mean ±2 SD, SD= 0.0001, 0.00003) and analytic results (solid lines). Fig. 1a shows the spiking probability in the bins lacking the signaling input from the simulation study with three levels of scaled diffusion coefficient (D/⌧ m V 2 ✓ = 0.004, 0.02, 0.04 ) as a function of the strength of the signaling input (circles and squares with dashed lines). The solid lines show the analytical spontaneous spiking probability (spiking probability without signaling input). The spiking probability in the bins lacking the signaling input is decreased from the analytic spontaneous spiking probability of a neuron (solid line). This deviation is stronger when the signaling inputs' amplitude is larger. We also calculated the spiking probability after the signal arrival, F A ( ). We counted the number of spikes of a postsynaptic neuron in the bins where the signal arrived and divided it by the total number of bins during which the signal occurred. Supplementary Fig. 1b shows F A ( ) increases for some values of excitatory signaling input.
We simulated two and three LIF neurons in different motifs to compare the pairwise and triplewise interactions determined from the simulation with analytic results (Supplementary Fig. 2). For the pairwise and triple-wise interactions between two and three neurons for excitatory shared inputs ( Supplementary Fig.2, top row), the analytical solutions ( Fig. 2c and Fig. 4c) capture the increasing and saturating trend of interactions as a function of the signal amplitude. There are, however, some discrepancies between the analytical and simulation results. They are caused by the theoretical mixture model that does not consider the carry-over effect ( Supplementary Fig. 1). In the simulation, we expected the spiking probability in bins lacking the signal to decrease, and the spiking probability after the signal arrival, to increase due to the carry-over effect. This effect would lead to stronger pairwise and triple-wise interactions than in the analytic results for excitatory shared inputs. We performed the same analysis on the inhibitory signaling inputs, where we observed more significant discrepancies due to the carry-over effect ( Supplementary  Fig.2, bottom row). However, we should note that the magnitude of the interactions caused by the inhibitory signaling inputs is much smaller than that of the excitatory signaling inputs. Further, the analytical method overestimates the interactions. Therefore, our conclusions in this study are not affected by these discrepancies.

Supplementary Note 2: Why the information-geometric measure is suitable for quantifying pairwise and higher-order interactions
We analyzed the pairwise and higher-order (i.e., triple-wise) statistical dependency among two and three neurons by using an information-geometric measure obtained from the neurons' interactions which are expressed in the form of the exponential family distribution (Eqs. 14 and 16) 37,62,64 . Here, we explain why we employed this measure instead of the alternatives.
For two neurons (Eq. 14), the information-geometric measure ✓ 12 of the pairwise interaction expresses their dependency, similarly to the Pearson's correlation coefficients computed for the binary patterns. Although zero values in both measures indicate that the two neurons are uncorrelated, their mathematical expressions are different. The information-geometric measure of the pairwise interaction can take an unconstrained real value whereas the Pearson's correlation coefficients are bounded in 1 and 1.
It is known that information-geometric measures are statistically more suitable for extracting pairwise and higher-order dependencies of the binary representation of neural activity 62,64 , compared with conventional measures such as Pearson's correlations, pairwise and triple-wise cross-correlations 4 , cumulants 5 , and mutual information 58,62,64 . For example, in the two-neuron system, the estimation of ✓ 12 from the data is not affected by the estimation of the activity rates of individual neurons. This property does not hold for the other measures such as the classical Pearson's correlation coefficient 64 . Namely, the simultaneous estimation of the firing rates and the neurons' dependency measured by the Pearson's correlation coefficients are correlated. On the contrary, the information-geometric measure of the correlation ✓ 12 is independent of the estimated firing rates of neurons. In this sense, it quantifies the genuine interaction that we cannot obtain from the activities of individual neurons: it can be inferred only if we simultaneously observe the two neurons. In information geometry, this property is known as orthogonality of the interaction to the activity rates of neurons.
The interactions among the three neurons (triple-wise interaction, ✓ 123 ) are obtained from Eq. 16. Similar to the information-geometric measure of the pairwise interactions, these triplewise interactions estimated from the spike data are not correlated with the lower-order statistics, namely the activity rates of individual neurons and the joint activity rates of two neurons. It thus quantifies a genuine triple-wise activity of the three neurons that can be inferred only if we observe the three neurons simultaneously.
The meaning of the non-zero triple-wise interaction would become clear if we compare it with the neural activity that expresses no triple-wise interaction. Here, suppose the triple-wise interaction is fixed at zero. In that case, the model corresponds to the maximum entropy model, which is the most unstructured model when the individual firing rates and joint firing rates of pairs of neurons (equivalently, pairwise correlations) are given. Such a model produces joint activities or inactivity of all three neurons (i.e., pattern '111' or '000') with non-zero probability, but these outcomes occur as a chance coincidence given the neurons' firing rates and pairwise correlations. The non-zero information-geometric measure of the triple-wise interaction ✓ 123 precisely measures the deviation of the joint activity or inactivity of the three neurons from this expected level. If it is positive, the three neurons are jointly active more frequently than the chance level expected from their activity rates and pairwise correlations. If it is negative, the neurons are simultaneously silent more frequently than the chance level.
A comparison of the other methods, such as covariance and mutual information (or Kullback-Leibler (KL) divergence), with the information-geometric measure is given in Amari 64 . The covariance method depends on firing rates, while the KL-divergence cannot reveal the sign of the correlation. Ohiorhenuan and Victor compared the KL divergence and the information-geometric measure, and their results pointed to the superiority of the information-geometric measure in practice 58 . Information-geometric measure can distinguish whether the interactions are more or less synchronous compared with that of chance level, while the KL divergence can measure only the magnitude. The information-geometric measure is also easier to calculate the confidence bound, and it was reported that it is robust to spike-sorting errors 58 . In summary, although information-geometric measure requires the data to be binarized, its statistical and practical advantages make it suitable for measuring the dependency of a neural population.

Supplementary Note 3: High-amplitude boundary for pairwise and triple-wise interactions induced by strong common excitatory and inhibitory inputs
Here, we investigated the case of an extreme strong shared signaling input and how it affects the pairwise and triple-wise interactions. This analysis helped us to find the high-amplitude boundaries for motifs in Fig. 6a. The first boundary line in the ✓ 123 versus ✓ 12 interaction plane for each motif arises from the high-amplitude limit (i.e., F A ( ) ! 1 for common excitatory inputs and F A ( ) ! 0 for common inhibitory inputs), while the spontaneous spiking probability F 0 ( ) changes within [0, 1]. These boundaries are shown in Fig. 6a as solid gray lines. Below, we explain, for each motif, how we found the relations for high-amplitude boundaries. First we consider pairwise interactions; then we will extend the formalism to triple-wise interactions for common excitatory and inhibitory inputs in star and (symmetric/asymmetric) triangle architectures.

Supplementary Note 3-a: Pairwise interaction for high-amplitude common input
Combining Eq. 15, Eq. 1, and Eq. 2, the pairwise interaction reads: This equation indicates that the common input induces a positive interaction for both excitatory and inhibitory connections. If the postsynaptic neuron receives an extremely strong inhibitory signal, the chance of its spiking in a short time window, , immediately diminishes after signal arrival: F A ( ) ' 0. This reduces Eq. S.1 to where a = b/(1 b) = /(1 ). Conversely, for an extremely strong excitatory signal, the chance of spiking in the time window, , is almost 1 immediately after signal arrival; hence F A ( ) ' 1. The pairwise interaction then simplifies to: Comparing Eq. S.2 and Eq. S.3, we can see that 1 F 0 ( ) simply replaces F 0 ( ). In Eq. S.3, for excitatory signaling input and with a > 0, we obtain a large ✓ 12 when F 0 ( ) ⌧ 1. In fact, ✓ 12 indefinitely increases as F 0 ( ) ! 0 + . On the other hand, for strong inhibitory signals in Eq. S.2, this criteria changes to 1 F 0 ( ) ⌧ 1; which means we need F 0 ( ) . 1 to obtain a high ✓ 12 . Consequently, in the low-activity regime, i.e., F 0 ( ) ⌧ 1, strong excitatory signals produce large ✓ 12 ; whereas in the high-activity regime, i.e., F 0 ( ) . 1, strong inhibitory inputs produce large pairwise interactions. The only assumption in the above reasoning is the low firing rate of the presynaptic signaling input: ⌧ 1. More accurately, the Poissonian probability of one or more signals arriving in a time window, , as long as  1, is b = 1 exp( ), which simplifies to b ' in the low firing rate limit of the signaling input. Accordingly, the prefactor, a, is in general, a = {1 exp( )}/ exp( ) which simplifies to a ' /(1 ) in the low input-rate limit. However, this correction does not affect the overall result of our analysis. Supplementary Fig. 3a shows how the pairwise interaction depends on the spontaneous spiking rate and the firing rate of the signaling input in a small time window, for inhibitory (left panel) and excitatory inputs (right panel). As the firing rate of the signaling input increases ( we still have  1), the amount of interaction induced by both excitatory and inhibitory inputs increases as well. However, increasing the spontaneous spiking probability of F 0 ( ), causes the interaction induced by the excitatory-to-pairs motif to decrease while that of inhibitory-to-pairs motif to increase.
Supplementary Note 3-d: High-amplitude for very small pairwise and triplewise interactions of excitatory/inhibitory to trio motif In the case of giving high-amplitude excitatory inputs to three neurons (F A ( ) ! 1), Eq. S.3 and Eq. S.7 draw the high-amplitude boundary over the whole range of F 0 ( ) in the ✓ 123 versus ✓ 12 plane, but a small region (0 < ✓ 12 < log(1 + a), 0 < ✓ 123 < log(1 + a)) remains, where the neuron operates in very high spontaneous regime, F 0 ( ) ! 1. Since the excitatory input is in high-amplitude regime ((F A ( ) ! 1) and simultaneously the spontaneous rate is very high (F 0 ( ) ! 1), the relations for pairwise and triple-wise interactions (Eq. S.1 and Eq. S.4) are indeterminate. In this regime, we can take the limit of the pairwise and triple-wise interactions (Eq. S.1, Eq. S.4 and Eq. S.5), by considering F A ( ) ! 1 and F 0 ( ) ! 1 simultaneously.

b
Strong excitatory input to pairs Strong inhibitory input to trio Strong inhibitory common input The pairwise interactions of two neurons generated by strong common inhibitory (left) and excitatory (right) inputs as a function of the spontaneous spike rate of a postsynaptic neuron F 0 ( ) and input rate . Common Inhibitory inputs generate strong pairwise interactions if the spontaneous rate of the postsynaptic neuron is high (i.e., F 0 ( ) ' 1), while common excitatory inputs generate strong interactions in the low spontaneous rate regime, F 0 ( ) ⌧ 1. (b) The triple-wise interactions among three neurons generated by strong common inputs in the inhibitory-to-trio (left) and excitatory-to-pairs (right) motifs. The same tendency found in the pairwise interactions applies to the triple-wise interaction. The excitatory-to-pairs motif generates strong triple-wise interactions in the low spontaneous rate regime, but the inhibitory-to-trio motif cannot induce such strong interactions in this regime.
Supplementary Note 4: Defining the analytical low (high) spontaneous rate boundaries for the excitatory (inhibitory) to trios/pairs in ✓ 123 versus ✓ 12 plane Here we describe in brief how to determine the low (high) spontaneous rate boundary for excitatory (inhibitory) input based motifs in the ✓ 123 versus ✓ 12 plane. One can easily find out how the spontaneous postsynaptic neuron's firing rate relates to the level of background synaptic activities from a stationary distribution of voltage trajectories for the LIF model neuron 8,100 : i.e., (Appendix V in Shomali et al. 36 ) where ⇥(V ) is the Heaviside step function: ⇥(V ) = 1 for V > 0 and otherwise ⇥(V ) = 0. µ is the mean firing rate: µ = ( R 1 0 sJ 0 (s) ds) 1 and V r is the resting voltage that we assumed to be zero.Ī is the mean current which is equal to V ✓ for the threshold regime. Since the stationary distribution is normalized ( , µ can be obtained explicitly 8 : Using Eq. S.24, each level of background input noise (diffusion coefficient) gives us a specific spontaneous firing rate that we can use to calculate F 0 ( ). In the star architecture with common excitatory input, for fixed F 0 ( ), by changing F A ( ) within the valid range [F 0 ( ), 1], we determined the whole range of triple-wise and pairwise interactions for excitatory-to-trio motif (i.e., Eq. S.15 and Eq. S.4) associated with a specific spontaneous rate (green to purple dashed lines in Fig.6a). We observed that, as the background input noise decreased, the region for excitatory-to-trio motif expands, so we used the lowest plausible spontaneous rate to determine the boundary (we used µ = 1 Hz for D = 5.5 ⇥ 10 17 ms · mV 2 ). This gave us the spontaneous spiking probability, F 0 ( ) = 0.005 (for = 5 ms). By changing F A ( ) within the range [F 0 ( ), 1] in Eq. S.15 and Eq. S.4, we obtained the second boundary for the excitatory-to-trio motif (Fig.6a, purple line).
For inhibitory-to-trio motif, a higher spontaneous rate shows stronger interactions. Thus, the second boundary is the limit of the high spontaneous rate. Since we are restricted to a low spontaneous rate regime (i.e., F 0 ( )  0.5), we chose the background noise for F 0 ( ) = 0.5, that is (D = 190 ms · mV 2 and µ = 100 Hz for = 5 ms). Here, because of the inhibitory input, F A ( ) varies within the range [0, F 0 ( )] in Eq. S.15 and Eq. S.4, but this time for the high spontaneous rate (µ = 100 Hz, F 0 ( ) = 0.5). The two boundaries for the inhibitory-to-trio region (Fig.6a left, bottom) lie very close to each other so the region is narrow in the plane of triple-wise versus pairwise interactions for low spontaneous rate regime (i.e., 0  F 0 ( )  0.5).
For inhibitory-to-pairs motif, again the second boundary is the high spontaneous rate ( D = 190 ms · mV 2 and µ = 100 Hz, F 0 ( ) = 0.5 for = 5 ms). This boundary was found by simulating LIF neurons for D = 190ms · mV 2 , to find the F A ( ) and F 2A ( ) relation in the neuron model. Then, we calculated pairwise and triple-wise interactions by inserting the obtained relation between F A ( ) and F 2A ( )(Eq. 13, Eq. S.15 and Eq. 17).
Finally, the second boundary for the excitatory-to-pairs motif in the low spontaneous rate regime is obtained using very low diffusion (µ = 1 Hz for D = 5.5 ⇥ 10 17 ms · mV 2 ). As in the case of the inhibitory-to-pairs motif, the boundary depends on both F A ( ) and F 2A ( ) and consequently on the neuron model. However, the LIF simulation in this very low spontaneous rate does not give reliable results (D = 5.5 ⇥ 10 17 ms · mV 2 , ⌧ m = 10ms, V ✓ = 5mV ). Thus, we took a different approach; we derived the postsynaptic neuron's CDF after signal arrival (i.e., F A ( ) and consequently F 2A ( )) for the very low spontaneous rate regime: where P s (v) is the stationary distribution (Eq. S.22). We find the instantaneous spiking just after the signal arrives, and because of very low diffusion, the neuron doesn't generate spike for a while after signal arrival. When A ⌧m  V ✓ , the CDF is: For a specific spontaneous firing rate of the postsynaptic neuron, by using Eq. S.26 and Eq. S.27, pairwise and triple-wise interactions are calculated (in Eq. S.15 and Eq. 17) and drawn in Fig.6a. As shown in Figure 6a, the boundary in the very low spontaneous rate regime, is determined (µ = 1 Hz) for both symmetric and asymmetric excitatory-to-pairs motif(purple dashed lines in the excitatory-to-pairs region, Fig.6a). At high spontaneous rates (Fig. 10b), the boundary for the lowest spontaneous rate (µ = 100 Hz) of the excitatory-to-pairs motif, is found by simulating the LIF neuron model at a specific diffusion that relates to µ = 100 Hz, finding the F A ( ) and F 2A ( ) relation, and calculating interactions from Eq. S.15 and Eq. 17. So like one boundary of inhibitory-to-pairs motif, this boundary depends on the neuron model (here LIF) in Fig. 6a. Regarding these boundaries (Supplementary Note 3, 4), we build a guide map for inferring basic motifs from neural interactions. A guide map for neural interactions is presented for postsynaptic low spontaneous rate of µ = 1 100 Hz (Fig. 10a). For high spontaneous rate, µ = 100 1000 Hz, the guide map changes in favor of inhibitory inputs motifs for strong interactions (Fig. 10b). Supplementary Fig. 4 shows the excitatory inputs regions inducing weak interactions between postsynaptic neurons at high spontaneous rates. In this regime, the excitatory inputs' motifs induce small interactions, making their regions smaller than those found at low postsynaptic spontaneous rates (Fig. 6a). The only fixed parameters are the bin size = 5 ms and the presynaptic rate = 5 Hz.
Note that the carry-over effect (Supplementary Note 1) decreases the spontaneous rate after applying the signaling input. Therefore, the results of the pairwise and triple-wise interactions (Fig. 6a) shift to lines for a lower spontaneous rate.
The dependence on bin size is depicted in Supplementary Fig. 5; it shows how different bin sizes change the regions of the guide map for spontaneous rates of µ = 1 100 Hz. When larger bin sizes are chosen, the regions for the excitatory motifs in low spontaneous rate regime (µ = 1 100 Hz), shrink.
Supplementary Note 4-a: How does excitatory-to-trio motif produce negative values for ✓ 123 ?
We observed that when the spontaneous rate of the postsynaptic neuron is less than µ = 5 Hz, the triple-wise interaction in the motif of excitatory inputs given to three neurons, becomes negative for some values of common input's amplitude (Fig. 6a). We can rewrite the pattern probabilities and triple-wise interaction for excitatory inputs given to three neurons from Eq. S.4 and Eq. S.5 as Pattern probabilities for neurons with spontaneous rate of 1 Hz receiving common signaling input with rate = 5 Hz (excitatory-to-trio motif). If ' 3 2 is bigger than ' 1 , the triple-wise interaction becomes negative. In this condition, we have log ' 3 2 > log ' 1 for some ranges of amplitudes which leads to the negative triple-wise interactions (c, green dots mark the range), (d, e, f ) The pattern probabilities when the spontaneous rate is 5 Hz. At F 0 = 5 Hz, there is no region in which log ' 3 2 > log ' 1 that induces negative triple-wise interactions (f, bottom right).
The triple-wise interaction is and For the excitatory-to-trio motif, x varies within the range 1  x  1/F 0 ( ). The spontaneous rate of µ = 1 Hz for a time window = 5 ms can generate negative triple-wise interactions for some range of common input amplitude ( Supplementary Fig. 6, top panels). However the spontaneous rate of µ = 5 Hz generates positive triple-wise interaction over the whole range of common input amplitude ( Supplementary Fig. 6, bottom panels). The term that makes the triple-wise interaction, negative at some values of common input amplitude for low rates is ' 2 : it is the fraction of pattern probability P (1, 1, 0) over P (1, 0, 0) ( Supplementary Fig. 6, b). When ' 3 2 exceeds ' 1 , Eq. S.29 gives a negative triple-wise interaction for excitatory-to-trio motif ( Supplementary Fig. 6, compare top with bottom panels).
By investigation, we found that when ⌘ < 0.25, where ⌘ = F 0 ( )(1 F 0 ( ))/a 2/3 , the triple-wise interaction becomes negative in the excitatory-to-trio motif. Under this condition, a negative triple-wise interaction occurs when the cumulative density function In addition to these synapses, we injected presynaptic inputs conveying our shared signaling inputs with specific frequency and amplitude. The shared inputs were simple exponential synapses with a decay time constant of 0.1 ms, and applied to a 0.7 section of soma. Supplementary Table 1 shows some parameters used in simulation of multicompartmental neuron model of pyramidal neuron in L5 (https://bbp.epfl.ch/nmc-portal, L5:TTPC1:cADpyr232:1).

Parameter
Description  Supplementary Fig 7. Examples of synaptic activation sites from distinct presynaptic neurons on the dendritic branches of the L5 pyramidal neuron. Each panels shows the morphology of the L5 pyramidal neuron (TTPC1:cADpyr232:1). The red and orange circles are the locations of the excitatory and inhibitory synapses, respectively. The last panel shows the L5 neuron with all synaptic inputs. All the panels are generated by the NEURON program and the data provided by the blue brain project (https://bbp.epfl.ch/nmc-portal).
input with a rate of 5 Hz for 300 s (Fig. 6b, circles). The amplitude of the signaling input was within the range 0.1 to 4. We saved the spike times of the Poissonian input and used it as shared input to run simulations on two other postsynaptic neurons. To reduce the computational cost, we performed one serial session for 900 s by repeating the common input's spike times three times. Then, we calculated the interactions of these three neurons when they received one shared input on top of other synapses, and plotted them in the ✓ 123 versus ✓ 12 plane (Fig. 6b). We repeated the above procedure five times for Poissonian shared input with a specific amplitude. As shown in Fig. 6b, as the amplitude of shared input increases, interactions get stronger, and all the data fall within the boundaries for the excitatory-to-trio motif. We also simulated interactions among three neurons under spontaneous conditions, i.e., when they all received the m-type synaptic inputs except for the one shared input. The interactions for the spontaneous activity are nearly zero, as expected (the black diamond in Fig. 6b).
We used the same approach on the excitatory-to-pairs motif by generating three Poissonian common inputs, where each postsynaptic neuron received two of them, in addition to the m-type synapses from different layers. We calculated the pairwise and triple-wise interactions and repeated the simulation five times for each amplitude (efficacy) of the shared signal within a range of 0.1 to 4 (Fig. 6b, squares). The simulated results are located near or in the predicted narrow excitatory-to-pairs region. For a large shared-input amplitude, the interactions reach the high-amplitude line in excitatory-to-pairs motif.
Next, we ran simulations on the inhibitory-to-trio and inhibitory-to-pairs motifs by using a multicompartmental neuron model of L5 in NEURON. The implementation was the same as the one for excitatory inputs, but this time, we used negative amplitudes as shared inputs and activated all other independent inputs (Poissonian with a rate of 10 Hz) the neuron received in the model of the blue brain project (L5:TTPC1:cADpyr232:1). Here, large amplitudes of inhibitory shared inputs lead to divergent results and the program encounters errors and terminates (in the blue brain NEURON program, the voltage after strong inhibitory shot goes to negative infinity), so we limited ourselves to two signaling input amplitudes (0.4 and 0.5). The result ( Supplementary  Fig. 8) shows that inhibitory inputs are very noisy and small, but the inhibitory-to-trio motif generates average negative triple-wise interactions (yellow circle) and the inhibitory-to-pairs motif generates average positive interactions, as predicted (yellow circle).
We also investigated the dependence of triple-wise and pairwise interactions on the simulation length. Supplementary Fig. 9 shows simulation results of pairwise (panel a) and triple-wise interactions (panel b) for the duration of 5s to 300s. They converged after about 150s (i.e. before 300s used in the above simulations). However, the duration required to obtain the desired accuracy in estimating neuronal interactions depends on the specific neuron models and synaptic inputs.

Supplementary Note 6: Comparing cross-correlation with interactions of log-linear model for distinguishing the motifs
In this section, we examine the cross-correlation measure for two and three neurons to see how well it can distinguish different motifs in the pairwise (CC 2 ) and triple-wise (CC 3 ) cross-correlation plane. The cross-correlation (the Pearson correlation coefficient) is defined by normalizing the covariance between two neurons (i.e. hx 1 x 2 i hx 1 ihx 2 i): is the activity of neuron 1 (2) and the hx 1 i (hx 2 i) is the expectation of activity x 1 (x 2 ). Although the standard cross-correlation is designed to measure a linear correlation of real-valued signals, here we applied it to binary data to examine its performance in discriminating population activity patterns (i.e., x 1 and x 2 are either 0 or 1). The cross-correlation (Eq. S.32) simplifies to: By representing the activity rates of neurons by probabilities of patterns, the cross-correlation can be written as 63,64 : P (1, 1) (P (1, 1) + P (1, 0))(P (1, 1) + P (0, 1)) p (P (1, 1) + P (1, 0))(P (0, 0) + P (0, 1))(P (1, 1) + P (0, 1))(P (0, 0) + P (1, 0)) , (S.34) where hx 1 x 2 i = P (1, 1), hx 1 i = P (1, 0) + P (1, 1) and hx 2 i = P (0, 1) + P (1, 1). The cross-correlation among three neurons is: where x 1 , x 2 and x 3 are the activities of neurons 1, 2 and 3, respectively. This equation can be rewritten as: Using pattern probabilities, the cross-correlation among three neurons becomes:   Color codes show each motif's region that is distinguishable from the rest except for the regions where the triple-wise correlation for excitatory-to-pairs and excitatory-to-trio motifs, diverges. The divergence occurs when the denominator of CC 3 (Eq. S.39) approaches zero. The excitatoryto-trio motif induces strong CC 2 and CC 3 as the amplitude of signaling input gets stronger. The inset is a magnification of the regions for the symmetric and asymmetric excitatory-to-pairs motifs. These motifs can induce strong pairwise but small triple-wise correlations in CC 3 versus CC 2 plane. (b) The inhibitory-to-trio motif induces a weak pairwise but strong triple-wise correlations. The spontaneous rate ranges within 1 100 Hz. Fixed parameters are = 5 ms and = 5 Hz.
Assuming identical neurons, the cross-correlation among three neurons can be simplified as: Using the cross-correlations CC 2 and CC 3 , we drew the boundaries that define the realm of correlations achievable under different motifs like the boundaries we defined before in the plane of the information-geometric parameters (i.e., ✓ 12 and ✓ 123 , Fig.6a). Supplementary Fig. 10 shows the boundaries in the CC 3 versus CC 2 plane for low spontaneous rate regime (i.e., 1 100 Hz); each motif has two boundaries: one arises from high-amplitude limit and the other from low (high) spontaneous rates for excitatory(inhibitory) inputs motifs. The motifs occupy distinct regions similar to the parameters of the log-linear model. In the cross-correlation plane, the excitatory-to-pairs motifs show strong pairwise cross-correlation but very weak triplet correlations even for large common signaling input amplitudes. In contract, inhibitory-to-trio motif shows stronger triplet but weaker pairwise correlation. Therefore, the regions are separable. However, we observed that in some range of F 0 t 0.5, the denominator of CC 3 goes to zero when s = 0.5 in Eq. S.39. Therefore CC 3 diverges near this point for high-amplitude boundaries of the excitatoryto-trio and excitatory-to-pairs motifs: this makes it difficult to identify the excitatory-based motifs when the spontaneous rates are near F 0 t 0.5. If the correlation from one measure is available, it is possible to go from the parameters of this measure to another measure. For example, cross-correlation can be written in terms of the parameters of the log-linear model. By inserting pattern probabilities from Eq. 14 in Eq. S.34, the cross-correlation as a function of interaction parameter, ✓ 12 becomes: If the two neurons are homogeneous (identical), the cross-correlation simplifies to: , (S.41) Supplementary Fig. 11 shows how the cross-correlation depends on the natural parameters of ✓ 1 and ✓ 12 (pairwise interaction). It is also possible to plot ✓ 12 as a function of the cross-correlation. The firing rates of neuron 1 and neuron 2 are µ 1 = P (1, 1) + P (1, 0) and µ 2 = P (1, 1) + P (0, 1).

Supplementary Note 7: How does adaptation modulate the model dependent boundary?
Here we examined more physiologically plausible extensions of the LIF neuron model to see how they affect the interaction parameters and boundaries in the triple-wise versus pairwise interaction plane. In particular, we added an adaptation effect to the leaky integrate-and-fire neuron model 14,72 : where the adaptation current is w, the coupling term of voltage to adaptation is a, and ⌧ w is the time constant of the adaptation variable. When the neuron generates a spike, the voltage is reset to V r and the adaptation variable w is increased by a parameter (b) for the spike-triggered adaptation, i.e., w = w + b. We used adaptation parameters (Supplementary Table 2) consistent with Brette and Gerstner 72 . Supplementary Fig. 13 shows the results with and without the adaptation at the threshold regime for the excitatory-to-pairs (rectangles) and inhibitory-to-trio (circles) motifs and three scaling amplitudes of signaling inputs, i.e., A/⌧ m V ✓ = 0.05, 0.1, 0.25. For adaptation, we considered two cases in which the coupling between voltage and adaptation current (a) is present or absent. The former (i.e. a > 0) considers subthreshold adaptation in addition to spike-triggered adaptation (b), while the latter (i.e. a = 0) assumes only the effect of spike-frequency adaptation.  Fig. 13 shows how adaptation modulates the pairwise and triple-wise interactions for inhibitory-to-trio (panel a) and excitatory-to-pairs motifs (panel b). Here, the adaptation makes the triple-wise (pairwise) interactions for the excitatory-to-pairs motif more negative (positive), while it makes the triple-wise (pairwise) interactions in the inhibitory-to-trio motif less negative (positive), for strong signaling inputs. This effect induces a larger difference in triple-wise and pairwise interactions between these two motifs.

Parameter
The above results show adding the effect of adaptation clearly differentiates the triple-wise interactions between the excitatory-to-pairs and inhibitory-to-trio motifs in a way that the excitatory-to-pairs motifs can induce stronger negative triple-wise and positive pairwise interactions, while inhibitory-to-trio motifs induce weaker triple-wise interactions. This phenomena can be mainly attributed to decreasing the spontaneous rate of postsynaptic neurons in the presence of adaptation (Supplementary Fig. 14). We observed that adaptation changes the firing rate of the postsynaptic neuron receiving common inputs on top of noisy background inputs from 22 Hz (excitatory-to-pairs motif) and 19 Hz (inhibitory-to-trio motif) to around 5 Hz and 3 Hz respectively (see also Supplementary Fig. 15, panel b and c). When there is no common input, adaptation decreases the firing rate of postsynaptic neurons from 22 Hz to around 3 Hz (compare number 1 blue and red in Supplementary Fig. 15, panel b). Therefore, adaptation, by reducing the spontaneous rate, increases the strength of interactions in excitatory-to-pairs motifs but decreases it in inhibitory-to-trio motifs ( Supplementary Fig. 15, panel b and c).
Other interesting effects of adaptation are that excitatory-to-pairs motifs can generate the observed CV (coefficient of variation) in V1 neurons 74 , and that the CV becomes high when the signaling input amplitude is increased (Supplementary Fig. 15a, solid red line). These results are consistent with experimental studies reporting high variability for neurons 16 . This effect of adaptation in increasing the variability has also been observed experimentally 73 . On the other hand, the inhibitory-to-trio motif cannot generate such a high CV (Supplementary Fig. 15a, solid blue line).
Finally, Supplementary Fig. 16 shows how adaptation modulates the model-dependent boundaries in the triple-wise versus pairwise interaction plane. Here, we used the adaptive exponential integrate-and-fire model (aEIF) 14,72 in which a term, slop ⇥ exp((V V ✓ )/slop), is added to Eq. S.47 72 . Supplementary Fig. 16a shows how the spiking probability of a neuron within 5 ms upon signal arrival (i.e. F A ( = 5)) varies as a function of the scaled amplitude of signaling input for three levels of the scaled diffusion (noisy inputs) and three neuron models (LIF, aLIF, aEIF). Each plot starts from the spontaneous spiking probability (µ ). By increasing the scaled amplitude of the signaling inputs, F A ( ) gradually increases and eventually reaches 1 at high scaled amplitudes. When the scaled diffusion is 2.5 ⇥ 10 10 or less, the probability of spiking within 5 ms doesn't change for aLIF (superposition of solid red and dashed light red lines) or aEIF (superposition of solid green and dashed light green lines), so the result for D = 2.5 ⇥ 10 10 can be used for defining the low spontaneous rate boundaries in the ✓ 123 versus ✓ 12 plane. Supplementary Fig. 16b shows the boundaries of low spontaneous rate regime for excitatory-to-pairs architectures for the adaptive LIF (aLIF, Eq. S.47, red lines) and the adaptive exponential IF (aEIF, blue lines) neuron models. Each line starts at zero and, as the amplitude of common signaling inputs increases, it reaches the high amplitude boundary's curve at specific point. The point it reaches the high amplitude boundary depends on the spontaneous rate: µ = 1 Hz for aLIF and 6.8 Hz for aEIF. The high-amplitude boundaries by changing the neuron model remain the same because they are independent of neuron model (gray lines). The motifs regions are still separable and can be distinguished for more biological plausible models like aLIF and aEIF. The boundaries of excitatoryto-pairs, symmetric and asymmetric for adaptive LIF (aLIF, red) and adaptive exponential LIF (aEIF, blue) neuron model. The results for the low spontaneous rate limit are from a simulation study of aLIF and aEIF neuron models (i.e., D/(⌧ m V 2 ✓ ) = 2.5 ⇥ 10 10 ) and define the low spontaneous rate boundary for excitatory-to-pairs motifs, while of the high-amplitude boundary is independent of the neuron model (gray line). The brown dashed line shows the low spontaneous rate boundary for the excitatory-to-trio motif. The common parameters for the three neuron models in panel A are the membrane time constant, ⌧ m = 10 ms, threshold voltage, and V ✓ = 20 mV. For adaptive neuron models, the adaptation time constant is ⌧ w = 144 ms, coupling parameter is a = 0.0133, and spike-triggered average is b = 2.4 mV. For aEIF, the slope factor is 2 mV 72 .

Supplementary Note 8: How do interactions change if we go away from the threshold regime?
Here, we investigated how the results change when we go beyond the threshold regime. We performed simulations when the postsynaptic neurons are in the subthreshold and suprathreshold regimes and compared the results with simulations at the threshold regime. Supplementary  Fig. 17 shows the interactions of the simulations in the triple-wise versus pairwise interaction plane. The circles and rectangles in orange display the results when the postsynaptic neurons operate at the threshold regime. For the same parameters, those in yellow and pink are obtained when the mean inputs to the postsynaptic neurons reach 2 and 10 percent below the threshold ( I/Ī = 0.02, 0.1). The circles show the interactions obtained for the inhibitory-to-trio motifs, while the rectangles are for the excitatory-to-pairs motifs. The results are for three scaled amplitudes of the signaling input at A/⌧ m V ✓ = 0.2 (no line inside the symbols), 0.4 (one line), 1 (two crossed lines). These results show that, as the mean input goes away from the threshold in the subthreshold regime, the interactions for the excitatory-to-pairs motifs get stronger while those for the inhibitory-to-trio get weaker.
These results can be interpreted as follows. At threshold regime, the mean input to a postsynaptic neuron is set very close to the threshold. Thus, even a small amount of noise induces spikes of the postsynaptic neuron. In subthreshold regime where the mean input given to the postsynaptic neuron is set below the threshold, the neuron generates a spike as long as the variance of the noise helps the voltage to reach the threshold. In this case, the postsynaptic neuron fires less frequently (smaller amount of F 0 ( )) compared with the case of being at the threshold regime. Since the neuron is in a regime of lower spontaneous rate, the interactions induced by the excitatory inputs increase while the interactions caused by inhibitory inputs decrease (see Supplementary Note 3). Therefore, the hypothesis that the excitatory-to-pairs motif is behind the empirically observed strong triple-wise and pairwise interactions is not only unchanged but also more strongly supported when the neuron operates in the subthreshold regime rather than operating at the threshold regime.
We then moved on to the suprathreshold regime, where the mean input is set above the threshold of the postsynaptic neuron's voltage ( Supplementary Fig. 17). The symbols marked in green and blue were obtained when the mean inputs to the postsynaptic neurons reached 2 and 10 percent above the threshold ( I/Ī = 0.02, 0.1). In the suprathreshold regime, the picture mentioned above reverses. Since the mean input is set above the threshold, a neuron generates regular spikes with high frequency, changing the regime of spontaneous spiking probability from low to high firing rates, F 0 ( ) > 0.5. In this regime, the inhibitory signaling input can induce stronger interactions while excitatory signaling input cannot. We reconfirmed this interpretation in the simulation study ( Supplementary Fig. 17). For the range of scaling amplitudes of signaling inputs A/⌧ m V ✓ = 0.2, 0.4, 1, we observed that the interactions induced by inhibitory inputs get stronger as the mean input increases to 2 and 10 percent above the threshold (compare the orange circles with the green and blue circles of the symbols with crossed lines, for example). The interactions for the motif of excitatory-to-pairs of neurons. When a neuron is in the subthreshold regime, the interactions are reduced for inhibitory-to-trio motif (compare the orange with the pink and yellow circles) while they are increased for those of excitatory-to-pairs motif (compare orange with pink and yellow rectangles). In the suprathreshold regime, the excitatory-to-pairs motif generates weaker interactions (green and blue rectangles), while the inhibitory-to-trio motif generates stronger interactions, especially for stronger input (green and blue circles with crossed lines inside). Fixed parameters are = 10 ms, D/⌧ m V 2 ✓ = 0.002, and = 5 Hz. Each symbol is the average of at least 10 10 steps of the run so as to keep the mean squared error on the order of 10 3 -10 4 .
Supplementary Note 9: Does mixing motifs, excluding the excitatory-to-pairs motif, generate strong negative triple-wise interactions?
We have four architectures, common inputs given to pairs or to trio of postsynaptic neurons with either common excitatory or inhibitory inputs. Based on our analysis in this paper on common inputs and also recurrent activity, we know that the excitatory-to-pairs motif, if it exists, can induce strong interactions of ✓ 12 > 0 and ✓ 123 < 0. However, we can ask the following question: by mixing the other motifs (i.e, inhibitory-to-trio, excitatory-to-trio and inhibitory-to-pairs), can we reach a large magnitude of ✓ 12 (> 0) and ✓ 123 (< 0)? For this purpose, we mixed the three motifs, two by two, to see whether such a case exists. We mixed two motifs with firing rates of 1 = 5 Hz and 2 = 5 Hz for a range of common input amplitudes in each motif (i.e. from A = 0 to A = 50). The results of mixing the excitatory-to-trio motif with the inhibitory-to-trio motif are shown in Supplementary Fig. 18a. Although the inhibitory-to-trio motif induces negative triple-wise interactions, the interactions are mostly positive when it is mixed with the excitatory-to-trio motif. The only region in which the mixing triple-wise interactions are negative is the weak or near-zero amplitudes of excitatory-to-trio motif ( Supplementary Fig. 18a). The mixture of excitatory-to-trio and inhibitory-to-pairs motifs induces positive triple-wise interactions ( Supplementary Fig. 18b) as expected from its each individual motif. The results for mixing the inhibitory-to-trio and inhibitory-to-pairs are shown in Supplementary Fig. 18c. Although the inhibitory-to-pairs motif induces positive triple-wise interactions, when it is mixed with the inhibitory-to-trio motif, negative triple-wise interactions occur for most amplitudes of mixing ( Supplementary Fig. 18c). The interactions induced in this mixing are in the range of interaction of the inhibitory-to-trio motif, which is not strong and does not interfere with the excitatory-to-pair motif's range of interaction. The triple-wise and pairwise interactions in each case, saturate at a specific amplitude of common input.  Inhibitory-to-trio and Inhibitory-to-pairs a b c Supplementary Fig 18. Effect of mixing the motifs of excitatory-to-trio, inhibitoryto-trio, and inhibitory-to-pairs. The results of interactions are shown for (a) the mixtures of excitatory-to-trio and inhibitory-to-trio, (b) excitatory-to-trio and inhibitory-to-pairs, and (c) inhibitory-to-trio with inhibitory-to-pairs (right). The amplitude of common input in each motif varied from 0 (i.e., the motif is inactive) to 50 (arrow shows increasing amplitude in each panel), and the interactions (✓ 123 and ✓ 12 ) were calculated respectively. Along each colored line (following the arrow in each panel), one motif's amplitude is fixed while the other mixed motif's amplitude varies from 0 to 50. For example, in panel a, for each colored line, the amplitude of the inhibitory-to-trio motif is fixed (see color box), while the amplitude of the excitatory-to-trio motif changes along the line from 0 to 50. The color box shows in each panel which motif has the fixed amplitude. Fixed parameters are = 10 ms, D/⌧ m V 2 ✓ = 0.003, ⌧ m = 10 ms, V ✓ = 5 mV, 1 = 5 Hz and 2 = 5 Hz.